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Abstract — Multiplicative noise (also known as speckle noise) 
models are central to the study of coherent imaging systems, 
such as synthetic aperture radar and sonar, and ultrasound and 
laser imaging. These models introduce two additional layers of 
difficulties with respect to the standard Gaussian additive noise 
scenario: (1) the noise is multiplied by (rather than added to) 
the original image; (2) the noise is not Gaussian, with Rayleigh 
and Gamma being commonly used densities. These two features 
of multiplicative noise models preclude the direct application of 
most state-of-the-art algorithms, which are designed for solving 
unconstrained optimization problems where the objective has 
two terms: a quadratic data term (log-likelihood), reflecting the 
additive and Gaussian nature of the noise, plus a convex (possibly 
nonsmooth) regularizer {e.g., a total variation or wavelet-based 
regularizer/prior). In this paper, we address these difficulties 
by: (1) converting the multiplicative model into an additive 
one by taking logarithms, as proposed by some other authors; 
(2) using variable splitting to obtain an equivalent constrained 
problem; and (3) dealing with this optimization problem using the 
augmented Lagrangian framework. A set of experiments shows 
that the proposed method, which we name MIDAL {multiplicative 
image denoising by augmented Lagrangian), yields state-of-the-art 
results both in terms of speed and denoising performance. 

Index Terms — Multiplicative noise, speckled images, total 
variation, variable splitting, augmented Lagrangian, synthetic 
aperture radar, Douglas-Rachford splitting. 



I. Introduction 

Many special purpose imaging systems use coherent de- 
modulation of reflected electromagnetic or acoustic waves; 
well known examples include ultrasound imaging, synthetic 
aperture radar (SAR) and sonar (SAS), and laser imaging. Due 
to the coherent nature of these image acquisition processes, the 
standard additive Gaussian noise model, so prevalent in image 
processing, is inadequate. Instead, multiplicative noise models, 
i.e., in which the noise field is multiplied by (not added to) 
the original image, provide an accurate description of these 
coherent imaging systems. In these models, the noise field is 
described by a non-Gaussian probability density function, with 
Rayleigh and Gamma being common models. 

In this introductory section, we begin by briefly recalling 
how coherent imaging acquisition processes lead to multi- 
plicative noise models; for a more detailed coverage of this 
topic, the reader is refeiTed to El, 123, O. We 

then review previous approaches for dealing with images 
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affected by multiplicative noise and finally briefly describe 
the approach proposed in this paper 

A. Coherent Imaging and Multiplicative Noise 

With respect to a given resolution cell of the imaging device, 
a coherent system acquires the so-called in-phase and quadra- 
ture components, which are the outputs of two demodulators 
with respect to, respectively, cos(wi) and sin(wt), where lo is 
the angular frequency of the carrier signal. Usually, these two 
components are collected into a complex number, with the in- 
phase and quadrature components corresponding to the real 
and imaginary parts, respectively [24]. The complex observa- 
tion from a given resolution cell results from the contributions 
of all the individual scatterers present in that cell, which 
interfere in a destructive or constructive manner, according to 
their spatial configuration. When this configuration is random, 
it yields random fluctuations of the complex observation, a 
phenomenon which is usually referred to as speckle. The 
statistical properties of speckle have been widely studied and 
are the topic of a large body of literature EH, 1^, |l29l, ||3l1 . 

Assuming a large number of randomly distributed scatterers 
and no strong specular reflectors in each resolution cell the 
complex observation is well modeled by a zero-mean complex 
Gaussian circular density {i.e., the real and imaginary parts 
are independent Gaussian variables with a common variance). 
Consequently, the magnitude of this complex observation fol- 
lows a Rayleigh distribution and the square of the magnitude 
(the so-called intensity) is exponentially distributed ||29l . ||3T| . 
The term multiplicative noise is clear from the following 
observation: an exponential random variable can be written 
as the product of its mean value, the so-called reflectance 
(the parameter of interest to be estimated) by an exponential 
variable of unit mean (the noise). 

The scenario just described, known as fully developed 
speckle, leads to observed intensity images with a character- 
istic granular appearance, due to the very low signal to noise 
ratio (SNR). Notice that the SNR, defined as the ratio between 
the squared intensity mean and the intensity variance, is equal 
to one (0 dB); this is a consequence of the equality between the 
mean and the standard deviation of an exponential distribution. 

B. Improving the SNR: Multi-look Acquisition 

A common approach to improving the SNR in coherent 
imaging consists in averaging independent observations of 
the same resolution cell (pixel). In SAR/SAS systems, this 
procedure is called multi-look (M-look, in the case of M 
looks), and each independent observation may be obtained by 
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a different segment of the sensor array. For fully developed 
speckle, the resulting averages are Gamma distributed and the 
SNR of an AI -look image is improved to M. 

Another way to obtain an M-look image is to apply a low 
pass spatial filter (with a moving average kernel with support 
size AI) to a 1-look fully developed speckle image, making 
evident the tradeoff between SNR and spatial resolution. This 
type of Af-look image can be understood as an estimate 
of the underlying reflectance, under the assumption that this 
reflectance is constant in the set of cells (pixels) included in the 
averaging process. A great deal of research has been devoted to 
developing space variant filters which average large numbers 
of pixels in homogeneous regions, yet avoid smoothing across 
reflectance discontinuities, in order to preserve details/edges 
||20 |. Many other speckle reduction techniques have been 
proposed; see |[3T| for a comprehensive review of the literature 
up to 1998. 

C. Estimation of Reflectance: Variational Approaches 

What is usually referred to as multiplicative noise removal 
is of course nothing but the estimation of the reflectance 
of the underlying scene. This is an inverse problem calling 
for regularization, which usually consists in assuming that 
the underlying reflectance image is piecewise smooth. In 
image denoising under multiplicative noise, this assumption 
has been formalized, in a Bayesian estimation framework, 
using Markov random field priors IT), ||3T1 . More recently, 
variational approaches using total variation (TV) regulariza- 
tion were proposed ffl, I!?), ijMl. Il37l, |[38l. In a way, these 
approaches extend the spatial filtering methods referred to in 
the previous subsection; instead of explicitly piecewise flat 
estimates, these approaches yield piecewise smooth estimates 
adapted to the structure of the underlying reflectance. 

Both the variational and the Bayesian maximum a posteriori 
(MAP) formulations to image denoising (under multiplicative, 
Gaussian, or other noise models) lead to optimization problems 
with two terms: a data fidelity term (log-likelihood) and a 
regularizer (log-prior). Whereas in Gaussian noise models, 
the data fidelity term is quadratic, thus quite benign from 
an optimization point of view, the same is no longer true 
under multiplicative observations. In 1 1 ], the data fidelity term 
is the negative log-likelihood resulting directly from the M- 
look multiplicative model, which, being non-convex, raises 
difficulties from an optimization point of view. Another class 
of approaches, which is the one adopted in this paper, also 
uses the AI -look multiplicative model, but yields convex data 
fidelity terms by formulating the problem with respect to the 
logarithm of the reflectance; see (B], ||27l, Jll], 1371, lEl, and 
references therein. A detailed analysis of several data fidelity 
terms for the multiplicative noise model can be found in ll38l . 

The combination of TV regularization with the log- 
likelihood resulting from the multiplicative observation model 
leads to an objective function, with a non-quadratic term (the 
log-likelihood) plus a non-smooth term (the TV regularizer), to 
which some research work has been recently devoted [I], ||27l . 
||34l . ||37ll . ||38]| . Even when the log-likelihood is convex (as in 
lIZTl . Il34l . IIJTI . If38l ). it does not have a Lipschitz-continuous 



gradient, which is a necessary condition for the applicability 
(with guaranteed convergence) of algorithms of the forward- 
backward spHtting (FBS) class H, in), gll. Methods 
based on the Douglas-Rachford splitting (DRS), which do not 
require the objective function to have a Lipschitz-continuous 
gradient, have been recently proposed l,12J . I.38J . 

D. Proposed Approach 

In this paper we address the (unconstrained) convex opti- 
mization problem which results from the AI -look multiplica- 
tive model formulated with respect to the logarithm of the 
reflectance. As shown in ll38l . this is the most adequate for- 
mulation to address reflectance estimation under multiplicative 
noise with TV regularization. We propose an optimization 
algorithm with the following building blocks: 

• the original unconstrained optimization problem is first 
transformed into an equivalent constrained problem, via 
a variable splitting procedure; 

• this constrained problem is then addressed using an 
augmented Lagrangian method. 

This paper is an elaboration of our previous work [5 |, where 
we have addressed multiplicative noise removal also via a 
variable splitting procedure. In that paper, the constrained 
optimization problem was dealt with by exploiting the recent 
split-Bregman approach ||23l , but using a splitting strategy 
which is quite different from the one in f2S\. 

It happens that the Bregman iterative methods recently 
proposed to handle imaging inverse problems are equivalent 
to augmented Lagrangian (AL) methods |[30l , as shown in 
|j42l. We prefer the AL perspective, rather than the Bregman 
iterative view, as it is a standard and more elementary opti- 
mization tool (covered in most textbooks on optimization). In 
particular, we solve the constrained problem resulting from the 
variable splitting using an algorithm (of the AL family) known 
as alternating direction method of multipliers (ADMM) II16II , 

EH, nza. 

Other authors have recently addressed the variational 
restoration of speckled images IfTI. IfTH, ||27l. Il3l. |[37l, Il38l. 
The commonalities and differences between our approach and 
the approaches followed by other authors will be discussed 
after the detailed description of our method, since this 
discussion requires notation and concepts which will be 
introduced in the next section. 

The paper is organized as follows. Section formulates 
the problem, including the detailed description of the mul- 
tiplicative noise model and the TV regularization we adopt 
to estimate the reflectance image. Section |III] reviews the 
variable splitting and the augmented Lagrangian optimization 
methods, which are the basic tools with which our approach 
is built. Section |IV] introduces the proposed algorithm by 
direct application of the basic tools introduced in Section |lll] 
Section[V]discusses related work. Section[Vl]reports the results 
of a series of experiments aiming at comparing the proposed 
algorithm with previous state-of-the-art competitors. Finally, 
Section IVIII concludes the paper 
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II. Problem Formulation 

Let y G R" denote an n-pixel observed image, assumed 
to be a sample of a random image Y, the mean of which 
is the underlying (unknown) reflectance image x e M", 
that is, E[Y] = x. Adopting a conditionally independent 
multiplicative noise model, we have 



Yi = X, N,, for i = 1, 



..,n. 



(1) 



where N G K" is an image of independent and identically 
distributed (iid) random variables with unit mean, E(A^i) = 1, 
following a common probability density function p^. In the 
case of Af-look fully developed speckle noise, is a Gamma 
density 

PN{n) 



(2) 



r(Af) 

thus with expected value M[N] = 1 and variance 

al^E[{N-E[N])'] 

Accordingly, we define the signal-to-noise ratio (SNR) asso 
dated to a random variable li, for i = 1, . . . , n, as 



SNR 



M. 



(3) 



An additive noise model is obtained by taking logarithms 
of O [T5l, ET), m, 1221, ESl. For an ai'bitrary pixel of the 
image (thus dropping the pixel subscript for simplicity), the 
observation model becomes 



log Y ^ log X + log N . 
G z w 

The density of the random variable W = log N is 

PwH = PNie^e"' 



(4) 



r(M) 



Mw - 

e e 



"M 



thus 



PG\zig\z) =pw{g- z). 



(5) 



(6) 



Invoking the conditional independence assumption, we are 
finally lead to 



l0gPG|z(g|z) 



\ogpw{gs - Zs) 



(7) 



= C-A/^(z, + e3=-^0, (8) 

s=l 

where C is a constant (irrelevant for estimation purposes). 

Using the MAP criterion (which is equivalent to a regular- 
ization method), the original image is inferred by solving an 
unconstrained minimization problem with the form 

z e argmini(z), (9) 

z 

where i(z) is the objective function given by 

L(z) = -logpG|z(g|z) + A<^(z). (10) 

n 

= (z, +eS=-^=) + A0(z) + A; (11) 



in (fTTI) . A is an irrelevant additive constant, cj) is the negative 
of the log-prior (the regularizer), and A is the so-called 
regularization parameter 

In this work, we adopt the standard isotropic discrete TV 
regularizer fSl, that is. 



z)^^ J(AJz)2 + (A^z)2, 



(12) 



where A^z and A^'z denote the horizontal and vertical first 
order differences at pixel s G {1, . . . , n}, respectively. 

Each term (zs + e^=~^=) of (fTTl i. corresponding to the 
negative log-likelihood, is strictly convex and coercive, thus so 
is their sum. Since the TV regularizer is also convex (though 
not strictly so), the objective function L possesses a unique 
minimizer [13] , which is a fundamental property, in terms of 
optimization. In contrast, the formulation of the problem in 
terms of the original variables x (rather than their logarithm) 
leads to a non-convex optimization problem yj, [.34J . As seen 
in [l], the uniqueness of the minimizer of that non-convex 
objective is not guaranteed in general. 

In this paper, we will address problem (|9]l using variable 
splitting and augmented Lagrangian optimization. In the next 
section, we briefly review these tools, before presenting our 
approach in detail. 

III. Basic Tools 

A. Variable Splitting 

Consider an unconstrained optimization problem in which 
the objective function is the sum of two functions, one of 
which is written as a composition: 



min /i(u) 

u6K" 



/2 (5(u)) . 



(13) 



mm 

u,veR" 

s.t. 



(14) 



Variable splitting is a very simple procedure that consists 
in creating a new variable, say v, to serve as the argument of 
/2, under the constraint that .g(u) = v. The idea is to consider 
the constrained problem 

/l(u)+/2(v) 

3(u) = V, 

which is clearly equivalent to unconstrained problem ( fT3] ). 
Notice that in the feasible set {(u, v) : ^(u) = v}, the 
objective function in ( fT4] i coincides with that in ( fT3] l. 

The rationale behind variable splitting methods is that it 
may be easier to solve the constrained problem (fl4] i than it 
is to solve its unconstrained counterpart (fT3] l. This idea has 
been recently used in several image processing applications 

0, (23], 1221, Eg. 

A variable splitting method was used in |40| to obtain a fast 
algorithm for TV-based image restoration. Variable splitting 
was also used in |4| to handle problems involving compound 
regularizers; i.e., where instead of a single regularizer \(f) 
as in ( fTTT l. one has a linear combination of two (or more) 
regularizers Ai0i+A2(^2- In HI, Il27l . and ||40l , the constrained 
problem (fl4] i is attacked by a quadratic penalty approach, i.e., 
by solving 

min /i(u) + /2(v) + ^ 



l.9(u) 



(15) 
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by alternating minimization with respect to u and v, while 
slowly taking /i to very large values (a continuation process), 
to force the solution of ( fTSl l to approach that of ( fT4b . which in 
turn is equivalent to ( fTSl l. The rationale behind these methods 
is that each step of this alternating minimization may be 
much easier than the original unconstrained problem ( fTSl ). The 
drawback is that as fi becomes very large, the intermediate 
minimization problems become increasingly ill-conditioned, 
thus causing numerical problems (see [30|, Chapter 17). 

A similar variable spUtting approach underUes the recently 
proposed split-Bregman methods f23l; however, instead of 
using a quadratic penalty technique, those methods attack the 
constrained problem directly using a Bregman iterative algo- 
rithm |42|. It has been shown that, when g is a linear function, 
i.e., the constraints in (fl4b are linear {e.g., g{u) — Gu), the 
Bregman iterative algorithm is equivalent to the augmented 
Lagrangian method |l42l, which is briefly reviewed in the 
following subsection. 

B. Augmented Lagrangian 

In this brief presentation of the augmented Lagrangian 
method, we closely follow |30|, which the reader should 
consult for more details. Consider an equality constrained 
optimization problem (which includes ( fT4l i as a particular 
instance, if g is linear) 

min -E(z) 

zGR<^ (16) 
s.t. Az b =0, 

where b e Rp and A e W""^, i.e., there are p linear equality 
constraints. The so-called augmented Lagrangian function for 
this problem is defined as 

£a(z, a, /i) = E{^) + A^(Az - b) + I ||Az - h\\l (17) 

where A G K^* is a vector of Lagrange multipliers and ijl>{) 
is the penalty parameter 

The so-called augmented Lagrangian method (ALM), also 
known as the method of multipliers ll26l . Il33l . consists in 
minimizing £^(z, A,/i) with respect to z, keeping A fixed, 
and then updating A. 

Algorithm ALM 

1. Set fc = 0, choose /i > 0, and Ao. 

2. repeat 

3. Zfc+i G argmiiiz -Ca(z, Afc, ^) 

4. Afc+i ^ Afe + ^(Azfc+i - b) 

5. fc fc + 1 

6. until stopping criterion is satisfied. 

Although it is possible (even recommended) to update the 
value of ji in each iteration lU, ll30l . we will not consider 
that option in this paper. Importantly, unlike in the quadratic 
penalty method, it is not necessary to take /i to infinity to 
guarantee that the ALM converges to the solution of the 
constrained problem (fTSl l. 

Notice that (after a straightforward complete-the-squares 
procedure) the terms added to i?(z) in the definition of the 
augmented Lagrangian >C^(z, A, /i) in (fTTI i can be written as 



a single quadratic term, leading to the following alternative 
form for the ALM algorithm: 

Algorithm ALM ( version II) 

1. Set fc = 0, choose /i > 0, zq, and do. 

2. repeat 

3. Zfc+i Gargmin^£;(z) + f||Az-dfcI|^ 

4. dk+i ^ dfc - (Azfc+i - h) 

5. fc fc + 1 

6. until stopping criterion is satisfied. 

This form of the ALM algorithm makes clear its equivalence 
with the Bregman iterative method (see ['421). 

It has been shown that, with adequate initializations, the 
ALM algorithm generates the same sequence as a proximal 
point algorithm (PEA) applied to the Lagrange dual of problem 
( fT6b : for further details, see ll28l . 1351, and references therein. 
Moreover, the sequence {d^} converges to a solution of this 
dual problem and that all cluster points of the sequence {z^} 
are solutions of the (primal) problem ( fTSl l l28l . 

C. Augmented Lagrangian for Variable Splitting 

We now show how ALM can be used to address problem 
(fT4l) . in the particular case where g(u) = u (i.e., g is the 
identity function). This problem can be written in the form 
( fT6l ) using the following definitions: 

z= " , b = 0, A= [I -I], (18) 

and 

S(z) = /i(u) + /2(v). (19) 

With these definitions in place, steps 3 and 4 of the ALM 
(version II) can be written as follows 

I'^+A G argmin|/i(u) + /2(v) + ^||u-v-dfe||2| (20) 
Vfe+iJ u,v L 2 J 

dfc+i 4- dfc - (ufc+i - Vfe+i). (21) 

The minimization (f20b is not trivial since, in general, it 
involves non-separable quadratic and possibly non-smooth 
terms. A natural solution is to use a non-linear block-Gauss- 
Seidel (NLBGS) technique, in which ( f20| i is solved by alternat- 
ing minimization with respect to u and v. Of course this raises 
several questions: for a given dj., how much computational 
effort should be spent in approximating the solution of (f20li? 
Does this NLBGS procedure converge? Taking just one step 
of this NLBGS scheme in each iteration of ALM leads to 
an algorithm known as the alternating direction method of 
multipliers (ADMM) HI, ED, El (see also El, JM], ifMl l: 

Algorithm ADMM 

1. Set fc = 0, choose /i > 0, Vq, and do. 

2. repeat 

3. Ufc+i G argminu/i(u) + i^Wu-Vk - dk\\l 

4. Vfe+i G argmiuv /2(v) + f ||ufc+i - v - dfe||| 

5. dfc+i ^ dfc - (Ufc+i - Vfc+i - b) 

6. fc 4- fc + 1 

7. until stopping criterion is satisfied. 
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For later reference, we now recall the theorem by Eckstein 
and Bertsekas 1 16] in which convergence of (a generalized ver- 
sion of) ADMM is shown. This theorem applies to problems 
of the form ( fTSl ) with ^(u) = Gu, i.e., 



min/i(u) + /2(Gu), 



uG 



(22) 



of which ( fT4l i is the constrained optimization reformulation. 

Theorem 1 (Eckstein-Bertsekas, [16]): Consider problem 
(22), where G e M''^" has Ml column rank and 
/i : M" ^ (-00,00] and /2 : R'' (-00,00] 
are closed, proper, convex functions. Consider arbitrary 
fj. > andvo,do £ M'^. Let {rjk > 0, fc = 0,1,...} and 
{vk > 0, fc = 0, 1, ...} be two sequences such that 



r/k < 00 and 



k=Q 



t^fc < 00. 

fe=0 



Consider three sequences {uk G M", k — 0,1,...}, {v^ G 
M'^, fc = 0, 1, ...}, and {dk £ M"*, fc = 0, 1, ...} that satisfy 



Vk 



> 

Vk > 
dfe+i = 



Uk+i - argmin/i(u) + ^||Gu-Vfc-dA;| 
u 2 

Vfe+i - argmin/2(v) + ^\\Guk+i-v-dk\\l 
dk - (Gufe+i - Vfc+i). 



Then, if (2^ has a solution, say u*, the sequence {uk} to u*. If 
(221) does not have a solution, then at least one of the sequences 
{vfc} or{dk} diverges. 

Notice that the ADMM algorithm defined above generates 
sequences {uk}, {vfc}, and {d^} which satisfy the conditions 
in Theorem [1] in a strict sense (i.e., with rjk = Vk — 0). One 
of the important consequences of this theorem is that it shows 
that it is not necessary to exactly solve the minimizations in 
lines 3 and 4 of ADMM; as long as the sequence of errors are 
absolutely summable, convergence is not compromised. 

The proof of Theorem[T|is based on the equivalence between 
ADMM and the so-called Douglas-Rachford splitting method 
(DRSM) applied to the dual of problem The DRSM was 
recently used for image recovery problems in |12|. For recent 
and comprehensive reviews of ALM, ADMM, DRSM, and 
their relationship with Bregman and split-Bregman methods, 
see 113, 1351. 



IV. Proposed Approach 

To address the optimization problem (|9]l using the tools 
reviewed in the previous section, we begin by rewriting it as 



(z, u) — argminL(z,u) 

z,u 

S.t. Z = U, 



with 



L(z,u) = A/V(z, + 



s=l 



+ A0(u). 



(23) 
(24) 

(25) 



Notice how the original variable (image) z was split into a 
pair of variables (z, u), which are decoupled in the objective 
function 



A* II n 1 2 

z — z 

2Af" " 



The approach followed in 11381 also considers a variable 
splitting, aiming at the application of the ADMM method. 
However, the splitting therein adopted is different from ours 
and, as shown below, leads to a more complicated algorithm 
with an additional ADMM inner loop. 

Applying the ADMM method to the constrained problem 
defined by (l23Tl-(l25]l leads to the proposed algorithm, which 
we call multiplicative image denoising by augmented La- 
grangian (MIDAL). Obviously, the estimate of the image x is 
computed as x = e^'', component- wise. 

Algorithm MIDAL 

1. Choose ^ > 0, A > 0, Uo, and do. Set fc ^ 0. 

2. repeat 

3. z' ^ Uk + dk 

n 

4. Zfe+i ^ argmin^ [zs + e^^"''^) 

s=l 

5. u' ^ Zfe - dfe ^ ^ 

6. Ufc+i ^argmin-||u-u'f + -0(u). 

u 2 

7. dk+i ^ dk - {zk+i - Uk+i) 

8. fc ^ fc + 1 

9. until a stopping criterion is satisfied. 

The minimization with respect to z (line 4) is in fact a 
set of n decoupled scalar convex minimizations. Each of 
these minimizations has closed form solution in terms of 
the Lambert W function |14). However, as in we apply 
the Newton method, as it yields a faster (and very accurate) 
solution by running just a few iterations. 

The minimization with respect to u (line 6) corresponds 
to solving a £2— TV denoising problem with observed image 
u' and regularization parameter A//i or, equivalently, to com- 
puting the so-called Moreau proximity operator (see llT3l ) of 
{X/fi)(j), denoted prox^^^^^^ at u'; i.e., for 7 > 0, 

prox^^(u') = argmmi||u - u'|p + j(l){u). (26) 

We use Chambolle's algorithm IS) to compute prox^^, al- 
though faster algorithms could be applied pOl. As stated 
in Theorem [T] this computation does not have to be solved 
exactly as long as the Euclidian norm of the errors are 
summable along the ADMM iterations (and thus along the 
MIDAL iterations). 

Still invoking Theorem [T] and assuming that the sequences 
of optimization errors with respect to z (line 4 of MIDAL 
pseudo-code) and u (line 6 of MIDAL pseudo-code) are 
absolutely summable, then MIDAL convergence is guaranteed, 
because /i and /2 are closed proper convex functions, and 
G = I has full column rank. 

V. Comments on Related Work 

We will now make a few remarks on related work. Notice 
how the variable splitting approach followed by the ADMM 
method allowed converting a difficult problem involving a non- 
quadratic term and a TV regularizer into a sequence of two 
simpler problems: a decoupled minimization problem and a 
TV denoising problem with a quadratic data term. In contrast. 
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the variable splitting adopted in ||38l leads to an intermediate 
optimization that is neither separable nor quadratic, which is 
dealt with by an inner DRS iterative technique. 

TV-based image restoration under multiplicative noise was 
recently addressed in fTl]. The authors apply an inverse scale 
space flow, which converges to the solution of the constrained 
problem of minimizing TV(z) under an equality constraint on 
the log-likelihood; this requires a carefully chosen stopping 
criterion, because the solution of this constrained problem is 
not a good estimate. Moreover, as evidenced in the experi- 
ments reported in flS'l, the method proposed in fTl] has a 
performance far from the state-of-the-art. 

In II27I . a variable spUtting is also used to obtain an objective 
function with the form 

£'(z,u) = L(z,u) + a||z - u||^; (27) 

this is the so-called splitting-and-penalty method. Notice that 
the minimizers of E{z, u) converge to those of (|23]|-(|24]| only 
when a approaches infinity. However, since E{z, u) becomes 
severely ill-conditioned when a is very large, causing numer- 
ical difficulties, it is only practical to minimize E{z, u) with 
moderate values of a; consequently, the solutions obtained 
are not minima of the regularized negative log-likelihood ( fTTT i. 
Nevertheless, the method exhibits good performance, although 
not as good as the method herein proposed, as shown in the 
experiments reported below. 

As mentioned in Subsection II-CI in the approach followed 
in ||T], the objective function is con-convex. In addition to a 
lack of guarantee of uniqueness of the minimizer, this feature 
raises difficulties from an optimization point of view. Namely, 
as confirmed experimentally in lIZTl . the obtained estimate 
depends critically on the initialization. 

Finally, we should mention that the algorithmic approach 
herein pursued can also be interpreted from a Douglas- 
Rachford spHtting perspective HIl, lHH, OSj. In ||T2J, that 
approach was applied to several image restoration problems 
with non-Gaussian noise, including a multiplicative noise case, 
but not with the Gamma distribution herein considered. 

VI. Experiments 

In this section, we report experimental results comparing the 
performance of the proposed approach with those of the recent 
state-of-the-art methods introduced in [15J and [,27 J . We chose 
to focus on those two references for two reasons: (a) they 
report quantitative results and the corresponding implementa- 
tions are available; (b) experimental results reported in those 
papers show that the methods therein proposed outperform 
other recent techniques, namely the above mentioned |IT| and 
ll37i . as well as the (non-iterative) block-Stein thresholding of 

m. 

The proposed algorithm is implemented in MATLAB 7.5 
and all the tests were carried out on a PC with a 3.0GHz Intel 
Core2Extreme CPU and 4Gb of RAM. AU the experiments 
use synthetic data, in the sense that the observed image is 
generated according to ([T]l-(|2]l, where x is the original image. 
In Table H] we list the details of the 16 experimental setups 
considered: the original images, their sizes, and the maximum 



TABLE I 
Experimental setup 



Exp. 


Image 


Size 


M 


^rain 


2:max 


A 


1 


Cameraman 


256 X 256 


3 


0.03 


0.9 


4 


2 


Cameraman 


256 X 256 


13 


0.03 


0.9 


6.5 


3 


Lena 


256 X 256 


5 


0.03 


0.9 


4.5 


4 


Lena 


256 X 256 


33 


0.03 


0.9 


9.8 


5 


Siml 


128 X 128 


2 


0.05 


0.46 


5.5 


6 


Sim2 


300 X 300 


1 


0.13 


0.39 


3.5 


7 


Sim3 


300 X 300 


2 


0.03 


0.9 


3 


8 


Fields 


512 X 512 


1 


101 


727 


3.5 


9 


Fields 


512 X 512 


4 


101 


727 


4.5 


10 


Fields 


512 X 512 


10 


101 


727 


6.7 


11 


Mmes 


512 X 512 


1 


10-* 


255 


2 


12 


Mmes 


512 X 512 


4 


10-* 


255 


2.7 


13 


NImes 


512 X 512 


10 


10-* 


255 


4 


14 


Cameraman 


256 X 256 


1 


7 


253 


2.7 


15 


Cameraman 


256 X 256 


4 


7 


253 


4.5 


16 


Cameraman 


256 X 256 


10 


7 


253 


6.1 



and minimum pixel values (xmax and Xmin); the M values 
(which coincides with the SNR (O); the adopted value of A 
(the regularization parameter in (fTTI) ') for our algorithm. 

Experiments 1 — 7 reproduce the experimental setup used 
in 111 and ll27ll ; experiments 8-16 follow those reported in 
ifTSl . The 8 original images used are shown in Figure [T] The 
values of Xmax and x^i^ are as in (15] and [TtI, for comparison 
purposes. Notice the very low SNR values (M values) for most 
observed images, a usual scenario in applications involving 
multiplicative noise. 

The focus of this paper is mainly the speed of the algorithms 
to solve the optimization problem thus the automatic 
choice of the regularization parameter A is out of scope. 
Therefore, as in HI, ifTSll . and ll27l . we select A by searching 
for the value leading to the lowest mean squared error with 
respect to the true image. 

Assuming that conditions of Theorem [T] are met, MIDAL is 
guaranteed to converge for any value of the penalty parameter 
jj, > 0. This parameter has, however, a strong impact in 
the convergence speed. We have verified experimentally that 
setting fi — \ yields good results. For these reason, we have 
used this setting in all the experiments. 

MIDAL is initialized with the observed noisy image. The 
quality of the estimates is assessed using the relative error (as 
in 1271) 



and the mean absolute-deviation error (MAE) (as in ifTSl ) 

l|x-x||i 

MAE = 

n 

where x = cxp(z) and || • ||2 and || • ||i stand for the £2 and £1 
norms, respectively. As in |f27l|, we use the stopping criterion 

||xfc+i - Xfc||2 

llxfclb - ' 

with TO = 4 in experiments 1 to 7, as in 1,27 J , and to = 2 in 
experiments 8 to 16. 



TO APPEAR IN THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2010. 



7 




,xlO 



Objective function 



Fig. I. The 7 original images used in the 16 experiments: from top to bottom 
and left to right: Cameraman, Lena, Siml, Sim2, Sim3, Fields, and Nimes. 



A. Computing the TV Proximity Operator 

MIDAL requires, in each iteration, the computation of the 
TV proximity operator, prox..^^, given by ( l26b . for which we 
use Chambolle's fixed point iterative algorithm |8|. Aiming 
at faster convergence of Chambolle's algorithm, and con- 
sequently of MIDAL, we initialize each run with the dual 
variables (see 1 8 1 for details on the dual variables) computed 
in the previous run. The underlying rationale is that, as 
MIDAL proceeds, the images to which the proximity operator 
is applied get closer; thus, by initializing the computation of 
the next proximity operator with the internal variables of the 
previous iteration, the burn-in period is largely avoided. 



1.5 
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1 

0.5 
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10 

time (sec) 



15 
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Fig. 2. Evolution of the objective function jilt for the setting of Experiment 
1, using Chambolle's fixed point iterative algorithm |8| to compute the 
TV proximity operator with different initializations. Top plot: initialization 
with the dual variables computed in the previous iteration. Bottom plot: 
initialization with the dual variables set to zero. 



Another perspective to look at this procedure, already re- 
ferred to, is given by Theorem [T] which states that there is 
no need to exactly solve the minimizations in each iteration, 
but just to ensure the minimization errors along the iterations 
are absolutely summable. The fulfilment of this condition is 
easier to achieve with the proposed initialization than with 
a fixed initialization. Fig. |2] illustrates this aspect. For the 
setting of Experiment 1, it shows the evolution of the objective 
function ( fTTT i when the dual variables are initialized with the 
ones computed in the previous iteration (left hand) and when 
the dual variables are initialized to zero (right hand). All the 
curves in the top plot reach essentially the same value. Notice 
that MIDAL takes, approximately, the same time for a number 
of fixed point iterations between 5 and 20 to compute the TV 
proximity operator For a number of iterations higher than 20, 
MIDAL time increases because, in each iteration, it runs more 
fixed point iterations than necessary. In the plots in the right 
hand side, we see that the minimum of the objective function is 
never reached, although it can be approximated for large values 
of the fixed point iterations. Based on these observations, we 
set the number of fixed point iterations to 20 in all experiments 
of this section. 
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B. Results 

Table unreports the results of the 16 experiments. The times 
for our algorithm and that of |27 | are relative to the computer 
mentioned above. The numbers of iterations are also given, but 
just as side information since the computational complexity 
per iteration of each algorithm is different. The initialization 
of the algorithm of |27 | is either the observed image or the 
mean of the observed image; since the final values of Err are 
essentially the same for both initializations, we report the best 
of the two times. 

For experiments 8 to 16, we did not have access to the code 
of ifTSl . so we report the MAE and Err values presented in that 
paper According to the authors, their algorithm was run for a 
fixed number of iterations, thus the computation time depends 
only on the image size. The time values shown in TableHIlwere 
provided by the authors and were obtained on a MacBook Pro 
with a 2.53GHz Intel CoreDuo processor and 4Gb of RAM. 

In experiments 1 to 7, our method always achieves lower 
estimation errors than the method of |27|. Notice that the gain 
of MIDAL is larger for images with lower SNR, corresponding 
to the more difficult problems. Moreover, our algorithm is 
faster than that of 1271 in all the experiments, by a factor 
larger than 3. 

In all the experiments 8 to 16, our algorithm achieves lower 
MAE than the method of ifTSl . Concerning the relative error 
Err, our algorithm outperforms theirs in 5 out of 9 cases, there 
is a tie in two cases, and is outperformed (albeit by a very 
small margin) in two cases. 

Figures [3j SI ID and |6] show the noisy and restored images, 
for the Experiments 1 to 4, 5 to 7, 8 to 10, and 11 to 13, 
respectively. Finally, Fig. |7] plots the evolution of the objective 
function L{uk) and of the constraint function \\zk — Ufc||| 
along the iterations, for the Experiment 1 (Cameraman image 
and M — 3). notice the decrease of approximately 7 orders 
of magnitude of ||zfe — u^Hj along the 21 MIDAL iterations, 
showing that, for all practical purposes, the constraint ( l24b is 
satisfied. 

VII. Concluding Remarks 

We have proposed a new approach to solve the optimization 
problem resulting from variational (equivalently MAP) esti- 
mation of images observed under multiplicative noise models. 
Although the proposed formulation and algorithm can be used 
with other priors (namely, frame-based), here we have focused 
on total-variation regularization. Our approach is based on two 
building blocks: (1) the original unconstrained optimization 
problem was first transformed into an equivalent constrained 
problem, via a variable splitting procedure; (2) this constrained 
problem was then addressed using an augmented Lagrangian 
method, more specifically, the alternating direction method of 
multipliers (ADMM). We have shown that the conditions for 
the convergence of ADMM are satisfied. 

Multiplicative noise removal (equivalently reflectance esti- 
mation) was formulated with respect to the logarithm of the 
reflectance, as proposed by some other authors. As a conse- 
quence, the multiplicative noise was converted into additive 
noise yielding a strictly convex data term (i.e., negative of 
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Fig. 3. Left column: observed noisy images for Experiments 1 to 4 witli 
AI = 3, 13, 5, 33, respectively. Right column: image estimates. 



the log-likelihood function), which was not the case with 
the original multiphcative noise model. A consequence of 
this strict convexity, together with the convexity of the total 
variation regularizer, was that the solution of the variational 
problem (the denoised image) is unique and the resulting 
algorithm, termed MIDAL (multiplicative image denoising by 
augmented Lagrangian), is guaranteed to converge. 

MIDAL is very simple and, in the experiments herein 
reported, exhibited state-of-the-art estimation performance and 
speed. For example, compared with the hybrid method in ifTSll , 
which combines curvelet-based and total-variation regulariza- 
tion, MIDAL yields comparable or better results in all the 
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TABLE II 

Experimental RESULTS. Iter, Err, and MAE denote, respectively, the number of iterations, the relative error, and the mean 

ABSOLUTE-DEVIATION ERROR. THE TIMES ARE REPORTED IN SECONDS. THE TIME FOR flSl IS REFERRED TO A DIFFERENT MACHINE (SEE TEXT). 
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Exp 


Err 


MAe. 


Iter 


Time 


Err 
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Time 


Err 
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Time* 


1 


0.130 
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1 


10 


0.151 


1 13 


32 








2 


0.090 


0.025 
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g 


0.098 


1 15 


36 








3 


0.111 


U.U30 


11 

Zj 


10 


0.118 


133 


38 








4 


0.069 


n nil 


1 1 


9 


0.071 


182 


52 








5 


0.128 




lO 


1.4 


0.143 


165 


8 








6 


0.069 


0.012 


32 


21 


0.083 


166 


70 








7 


0.137 


0.023 


17 


11 


0.174 


165 


69 








8 


0.089 


29.06 


61 


121 








0.096 


32.67 


245 


9 


0.066 


20.86 


34 


68 








0.066 


22.0 


245 


10 


0.056 


17.94 


29 


58 








0.055 


18.24 


245 


11 


0.301 


12.38 


25 


14 








0.314 


13.27 


245 


12 


0.217 


8.78 


25 


50 








0.217 


8.98 


245 


13 


0.170 


6.87 


23 


46 








0.174 


7.11 


245 


14 


0.167 


12.74 


33 


15 








0.192 


16.78 


54 


15 


0.124 


9.43 


19 


9 








0.131 


10.67 


54 


16 


0.097 


7.42 


56 


26 








0.091 


7.44 


54 
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Fig. 6. Lett column: observed noisy images for Experiments 11 to 13. 
Right column: image estimates. Note: for better visualization, all the images 
underwent the nonUnear transformation (■)'''^ prior to being displayed. 



experiments. 

We are currently working on extending our approach to 
problems involving linear observation operators {e.g., blurs), 
other non-additive and non-Gaussian noise models, such as 
Poissonian observations [19J, [36] , and other regularizers, such 
as those based on frame representations. 
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